Method of determining subsidence in a reservoir

ABSTRACT

A reservoir simulator first estimates rock displacement parameters (u, v, and w) representing rock movement in the x, y, and z directions. When the rock displacement parameters (u, v, w) are determined, “ε x,y,z ” (the X,Y,Z elongation strains) and “γ xy,yz,zx ” (the shear strains) are determined since “ε x,y,z ” and “γ xy,yz,zx ” are function of “u”, “v”, and “w”. When “ε x,y,z ” and “γ xy,yz,zx ” are determined. “σ x,y,z ” (the elastic normal rock stress in x,y,z directions) and “τ xy,yz,xz ” (the elastic shear stress) are determined since “σ x,y,z ” and “τ xy,yz,xz ” are a function of “ε x,y,z ” and “γ xy,yz,zx ”. When “σ x,y,z ” and τ xy,yz,xz ” are determined, the rock momentum balance differential equations can be solved, since these equations are a function of “σ x,y,z ” and “τ xy,yz,xz ”. When any residuals are substantially equal to zero, the estimated rock displacement parameters (u, v, and w) will represent accurate rock displacement parameters for the reservoir. When the rock momentum balance differential equations are solved, the rock displacement parameters (u, v, w), at an advanced time, are known. These rock displacement parameters (u, v, w) represent and characterize a subsidence in a seabed floor because subsidence will result from rock movement; and rock movement will result from a continual withdrawal of oil or other hydrocarbon deposits or other fluids such as water over a period of time from an Earth formation. This abstract of the disclosure is given for the sole purpose of allowing a patent searcher to easily determine the content of the disclosure in this application.

BACKGROUND OF THE INVENTION

[0001] The subject matter of the present invention relates to areservoir simulator including a parameter determination software whichdetermines displacement parameters representing subsidence in anoilfield reservoir.

[0002] There are many recent reports of geomechanical modelling beingused predictively for evaluation of alternative reservoir developmentplans. In the South Belridge field, Kern County, Calif., Hansen et al¹calibrated finite-element models of depletion-induced reservoircompaction and surface subsidence with observed measurements. The stressmodel was then used predictively to develop strategies to minimiseadditional subsidence and fissuring as well as to reduce axialcompressive type casing damage. Berumen et al² developed an overallgeomechanical model of the Wilcox sands in the Arcabuz-Culebra field inthe Burgos Basin, northern Mexico. This model, combined with hydraulicfracture mapping together with fracture and reservoir engineeringstudies, was used to optimise fracture treatment designs and improve theplanning of well location and spacing.

[0003] The subject of fluid flow equations which are solved togetherwith rock force balance equations has been discussed extensively in theliterature. Kojic and Cheatham^(3,4) present a lucid treatment of thetheory of plasticity of porous media with fluid flow. Both the elasticand plastic deformation of the porous medium containing a moving fluidis analyzed as a motion of a solid-fluid mixture. Corapcioglu and Bear⁵present an early review of land subsidence modelling and then present amodel of land subsidence as a result of pumping from an artesianaquifer.

[0004] Demirdzic et al^(6,7) have advocated the use of finite volumemethods for numerical solution of the stress equations both in complexdomains as well as for thermo-elastic-plastic problems.

[0005] A coupling of a conventional stress-analysis code with a standardfinite difference reservoir simulator is outlined by Settari andWalters⁸. The term “partial coupling” is used because the rock stressand flow equations are solved separately for each time increment.Pressure and temperature changes as calculated in the reservoirsimulator are passed to the geomechanical simulator. Updated strains andstresses are passed to the reservoir simulator which then computesporosity and permeability. Issues such as sand production, subsidence,compaction that influence rock mass conservation are handled in thestress-analysis code. This method will solve the problem as rigorouslyas a fully coupled (simultaneous) solution if iterated to fullconvergence. An explicit coupling, i.e. a single iteration of the stressmodel, is advocated for computational efficiency.

[0006] The use of a finite element stress simulator with a coupled fluidflow option is discussed by Heffer et al⁹ and by Gutierrez and Lewis¹⁰.

[0007] Standard commercial reservoir simulators use a single scalarparameter, the pore compressibility, as discussed by Geertsma¹¹ toaccount for the pressure changes due to volumetric changes in the rock.These codes generally allow permeability to be modified as a function ofpore pressure through a table. This approach will not be adequate whenthe flow parameters exhibit a significant variation with rock stress.Holt¹² found that for a weak sandstone, permeability reduction was morepronounced under non-hydrostatic applied stress, compared with theslight decrease measured under hydrostatic loading. Rhett and Teufel¹³have shown a rapid decline in matrix permeability with increase ineffective stress. Ferfera et al¹⁴ worked with a 20% porosity sandstoneand found permeability reductions as high as 60% depending on therelative influence of the mean effective stress and the differentialstress. Teufel et al¹⁵ and Teufel and Rhett¹⁶ found, contrary to theassumption that permeability will decrease with reservoir compaction andporosity reduction, that shear failure had a beneficial influence onproduction through an increase in the fracture density.

SUMMARY OF THE INVENTION

[0008] Oil recovery operations are seeing increased use of integratedgeomechanical and reservoir engineering to help manage fields. Thistrend is partly a result of newer, more sophisticated measurements thatare demonstrating that variations in reservoir deliverability arerelated to interactions between changing fluid pressures, rock stressesand flow parameters such as permeability. Several recent studies, forexample, have used finite-element models of the rock stress tocomplement the standard reservoir simulation.

[0009] This specification discusses current work pertaining to: fullyand partially coupling geomechanical elastic/plastic rock stressequations to a commercial reservoir simulator. This finite differencesimulator has black-oil, compositional and thermal modes and all ofthese are available with the geomechanics option. In this specification,the implementation of the aforementioned stress equations into the code,hereinafter called the ‘Parameter Determination Software’, is discussed.Some work on benchmarking against an industry standard stress code isalso shown as well as an example of the coupled stress/fluid flow. Thegoal in developing this technology within the simulator is to provide astable, comprehensive geomechanical option that is practical forlarge-scale reservoir simulation.

[0010] Further scope of applicability of the present invention willbecome apparent from the detailed description presented hereinafter. Itshould be understood, however, that the detailed description and thespecific examples, while representing a preferred embodiment of thepresent invention, are given by way of illustration only, since variouschanges and modifications within the spirit and scope of the inventionwill become obvious to one skilled in the art from a reading of thefollowing detailed description.

BRIEF DESCRIPTION OF THE DRAWINGS

[0011] A full understanding of the present invention will be obtainedfrom the detailed description of the preferred embodiment presentedhereinbelow, and the accompanying drawings, which are given by way ofillustration only and are not intended to be limitative of the presentinvention, and wherein:

[0012]FIG. 1 illustrates control volume (light grey grid) for Rock ForceBalance;

[0013]FIG. 2 illustrates comparison of Rock Displacement Predictions(case 1);

[0014]FIGS. 3A and 3B illustrate comparison of Rock DisplacementPredictions (case 2);

[0015]FIGS. 4a 1, 4 a 2, 4 b 1, 4 b 2, 4 c 1, and 4 c 2 illustrate totalnormal stress in Mid X-Y plane;

[0016]FIG. 4d illustrates Table 1, and Table 1 further illustratessimulation parameters used in comparing ABAQUS and Reservoir SimulatorPredictions);

[0017] FIGS. 5-9 illustrate examples of structured and unstructuredgrids;

[0018]FIGS. 10 and 11a illustrate examples of a staggered grid;

[0019]FIGS. 11b, 12 and 13 illustrate the “Flogrid” software and the“Eclipse simulator software” and the generation of ‘simulation results’which are displayed on a 3-D viewer;

[0020]FIG. 14 illustrates an example of the ‘simulation results’;

[0021]FIG. 15 illustrates the fact that the “Eclipse simulator software”includes the ‘Parameter Determination Software’ of the presentinvention;

[0022]FIG. 16 illustrates the results which are generated when the‘Parameter Determination Software’ is executed by a processor of aworkstation, the results being presented and displayed on a recorder ordisplay device;

[0023]FIG. 17 illustrates how the results generated in FIG. 16 assist indetermining subsidence in an oilfield reservoir; and

[0024]FIG. 18 is a detailed flowchart of the ‘Parameter DeterminationSoftware’ of the present invention.

DESCRIPTION OF THE INVENTION

[0025] Elastic Stress Equations

[0026] Steady state “rock momentum balance equations” in the x, y and zdirections can be written as¹⁷: $\begin{matrix}{\begin{matrix}{{\frac{\partial\sigma_{x}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{x} + P_{x}} = 0} \\{{\frac{\partial\sigma_{y}}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{zy}}{\partial z} + F_{y} + P_{y}} = 0} \\{{\frac{\partial\sigma_{x}}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + F_{z} + P_{z}} = 0}\end{matrix}.} & (1)\end{matrix}$

[0027] Here {right arrow over (F)} is the body force, {right arrow over(P)} is the interaction force between the solid and the fluid. Thisinteraction force is discussed further below.

[0028] The elastic normal stresses σ and shear stresses τ can beexpressed in terms of strains, ε and γ as: $\begin{matrix}{\begin{matrix}{\sigma_{x} = {{2G\quad ɛ_{x}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\sigma_{y} = {{2G\quad ɛ_{y}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\sigma_{z} = {{2G\quad ɛ_{z}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\tau_{xy} = {G\quad \gamma_{xy}}} \\{\tau_{yz} = {G\quad \gamma_{yz}}} \\{\tau_{zx} = G_{\gamma_{zx}}}\end{matrix}.} & (2)\end{matrix}$

[0029] Constants G, also known as the modulus of rigidity, and λ areLamé's constants. They are functions of Young's modulus, E, andPoisson's ratio, v. Strains ε_(x,y,z) are defined in terms ofdisplacements in the x,y,z directions, namely u,v and w. Thus$\begin{matrix}{\begin{matrix}{ɛ_{x} = \frac{\partial u}{\partial x}} \\{ɛ_{y} = \frac{\partial v}{\partial y}} \\{ɛ_{z} = \frac{\partial w}{\partial z}} \\{\gamma_{xy} = {\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}}} \\{\gamma_{yz} = {\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}}} \\{\gamma_{zx} = {\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}}}\end{matrix}.} & (3)\end{matrix}$

[0030] Within the simulator, the rock force is calculated in units ofNewtons, lbf, or dynes depending on whether the user has chosen mks,field or cgs units for the simulation.

[0031] Gridding

[0032] Before returning to a discussion of the fluid component massconservation and rock mass conservation equations, the special griddingfor the rock force (momentum) balance equations is briefly outlined.

[0033] To implement these rock momentum balance equations in a finitedifference simulator, separate control volumes for balancing the forceare used in each orthogonal grid direction. These control volumes arestaggered in each of the coordinate directions.

[0034] In FIG. 1, an x-direction rock force balance control volume isshown. Also shown in this Figure is the finite difference discretizationused. The grid with the darker black lines in this figure is thestandard reservoir simulation grid which are control volumes for fluidcomponent and rock mass conservation. Rock displacements u, v and w aredefined on the faces of this grid. Rock momentum control volumes in aparticular coordinate direction are centred on the rock displacement inthat direction. In FIG. 1 is shown the x-momentum control volume (lightgrey grid) in the centre of which is the x-direction rock displacement,u, which is defined on the face between two control volumes for massconservation (black grid). Similar rock momentum control volumes are setup in the y and z directions around v and w, the y and z rockdisplacements

[0035] This grid is fixed in time. Rock will flow elastically orplastically in the mass conservation grid. A correction term to Darcy'slaw that accounts for the extra transport of fluid component mass intoan adjoining mass conservation cell due to rock movement into that cellcan be calculated in the simulator and is discussed below.

[0036] We have found the staggered gridding to be the most natural fordifferencing purposes and for setting of boundary conditions. It is alsothe most accurate because it achieves second order accuracy whencalculating the normal stresses and mixed first/second order accuracywhen calculating the shear stresses. To illustrate, if the user hasinput the more general stiffness matrix¹⁸, {right arrow over (S)}_(N),instead of the Young's moduli and Poisson ratios, then the differencednormal stress on the x+ side of the rock momentum balance control volumein FIG. 1 can be written $\begin{matrix}{\begin{matrix}{\sigma_{x +} = {{S_{N_{{xx}_{x +}}} \cdot ɛ_{x_{x +}}} + {S_{N_{xy}} \cdot ɛ_{y_{x +}}} + {S_{N_{xz}} \cdot ɛ_{z_{x +}}}}} \\{ɛ_{x_{x +}} = \frac{u_{{i + 1},j,k} - u_{i,j,k}}{\Delta \quad x_{i,j,k}}} \\{ɛ_{y_{x +}} = \frac{v_{i,{j + 1},k} - v_{i,j,k}}{\Delta \quad y_{i,j,k}}} \\{ɛ_{z_{x +}} = \frac{w_{i,j,{k + 1}} - w_{i,j,k}}{\Delta \quad z_{i,j,k}}}\end{matrix}.} & (4)\end{matrix}$

[0037] Another consequence of the staggered gridding is that the forcebalance equations (1) give rise to a thirteen-banded jacobian forthree-dimensional calculations. In addition to the standard diagonal andsix off-diagonal bands, there are six more connecting points in thedifferencing stencil which contribute the other bands.

[0038] Fluid Component Conservation, Rock Mass Conservation and VolumeBalance Equations

[0039] Compositional and thermal simulation require the user to specifythe number of fluid components. Fluid components are conserved in thestandard user-specified grid blocks as noted above. Rock also flowseither elastically or plastically, in addition to fluid, amongst thesegrid cells. A rock mass conservation equation can be written inpseudo-differenced form as: $\begin{matrix}{{V \cdot \rho_{rock} \cdot \left\lfloor {\left( {1 - \phi} \right)^{n + 1} - \left( {1 - \varphi} \right)^{n}} \right\rfloor} = {\sum\limits_{faces}{A_{face} \cdot \left( {u,v,{w_{face}^{n + 1} - u},v,w_{face}^{n}} \right) \cdot \left( {\rho_{rock} \cdot \left( {1 - \phi} \right)^{n + 1}} \right)^{up}}}} & (5)\end{matrix}$

[0040] where

[0041] V=bulk volume of the rock mass conservation cell

[0042] ρ_(rock)=specific density of rock

[0043] φ=porosity of rock (variable)

[0044] u,v,w=x,y,z total rock displacement

[0045] Superscripts n+1, n refer to the time level−n refers to the lasttime step, n+1 refers to the current or latest time level. Thesuperscript “up” refers to the upstream of the rock displacement changeover the time step. Flow of rock into an adjoining cell originates fromthe cell away from which rock flows during the timestep, as denoted byu,v,w_(face) ^(n+1)−u,v,w_(face) ^(n) in the above equation (5). Totalrock displacement is the sum of elastic and plastic rock displacement.The elastic rock displacements in the coordinate directions and porosityare included in the simulator variable set. Rock properties, for examplerelative permeability endpoints, are convected.

[0046] Fluid component mass conservation assumes the usual form and isdiscussed elsewhere^(19,20).

[0047] In order that the fluid volumes fit into the available porevolume, a final fluid phase volume balance equation is included. It isoften written as the sum of fluid phase saturations (volume of fluiddivided by pore volume) equal to unity.

[0048] Additional Enhancements to the Balance Equations

[0049] A term that describes the change of rock momentum in time, namely${m_{rock}\frac{\partial^{2}u}{\partial t^{2}}},$

[0050] where m_(rock) is the mass of rock in the rock balance grid cell,can also optionally be switched on by the user. By default it is omittedbecause this term is small and makes little contribution to the forcebalance. Also, when it is not specified, wave solutions of the formf(x±ct) are not calculated where c is some characteristic rockcompression/shear velocity. These wavelets may cause instability andusually are of little interest over standard simulation time scales.

[0051] Secondly, as discussed by Corapcioglu and Bear⁵, in a deformingporous medium it is the specific discharge relative to the moving solid,q_(r), that is expressed by Darcy's law

q _(r) =q−φV _(s) =−K·∇H   (6)

[0052] where

[0053] q=specific discharge of a fluid component (flux)

[0054] φ=porosity

[0055] V_(s)=velocity of the solid

[0056] K=fluid phase permeability

[0057] H=pressure head, including the gravitational head

[0058] This term is, again, omitted by default but can be included ifdesired. In problems where the elastic and plastic rock flow is smalland near steady state, there is not much contribution from this term. Insituations where the pressure field is changing continuously andstrongly, for example in a thermal cyclic huff and puff scenario wherethe near wellbore pressure can vary by more than 500 psi through a shortcycle time period, this term can have a strong effect on the simulation.In particular, it can reduce the unphysically high injection pressuresnormally predicted by a conventional simulator. There is an additionalCPU expense associated with this term since it can add one or morenewton iterations to a timestep.

[0059] Interaction Between the Fluid and the Rock

[0060] The interactive force between the fluid and solid, denoted{overscore (P)} in equation (1), can be expressed in the form⁵

{right arrow over (P)}={right arrow over (F)} _(f) ^(in) +φ{right arrowover (F)} _(f) +∇·{right arrow over (T)} _(f)   (7)

[0061] where

[0062] {overscore (F)}_(f) ^(in)=inertial force of the fluid per unitbulk volume

[0063] {overscore (F)}_(f)=body force acting on the fluid per unit fluidvolume

[0064] {overscore (T)}_(f)=stress tensor of the fluid

[0065] The stress tensor of the fluid can be related to the pore fluidpressure, following Terzaghi's²¹ concept of intergranular stress,

{right arrow over (T)} _(f) =−p·{overscore (I)}  (8)

[0066] where {right arrow over (I)} is the unit tensor. We neglect theeffect of the inertial forces of the fluid and write the interactionforce as

{right arrow over (P)}=−∇p+φ{right arrow over (F)} _(f)   (9)

[0067] The pore pressures are calculated at the centres of the massconservation cells so that the pore pressure gradient in equation (9),∇p, is centred in the rock momentum conservation control volume as itshould be.

[0068] Interaction of the rock stress on flow parameters, in particularfluid permeability, is a subject that is the topic of much currentresearch. As noted above in the introduction, many authors are nowcarrying out triaxial compression tests on core and rock samples inwhich fluid is flowing in the stressed rock matrix. Most commercialreservoir simulators contain permeability or transmissibility modifiersas a function of fluid pore pressure. This feature was used in ageomechanics simulation study by Davies and Davies²² where therelationship between stress and permeability was calculated based on therelationship between the stress and porosity, and between permeabilityand porosity for various rock types. Yale and Crawford²³ discuss thework of Holt¹², Teufel and Rhett^(13,15,16) and others. They also usedcritical state theory to model observed stress-permeability experimentaldata. It was concluded that this theory did well represent the observeddata.

[0069] We are currently exploring several constitutive relationshipsbetween stress and permeability. A mean effective stress versuspermeability is one of these relationships, another is the classicalpermeability-porosity. Another possibility is permeability versus σ₁−σ₃,the difference between the maximum and minimum principal stresses, asdisplayed on the web by GFZ-Potsdam²⁴. The success of the critical statetheory as discussed in the above paragraph is also of interest.

[0070] Variables and Solution of the Equations

[0071] We solve for moles of all fluid components per unit pore volume,energy per unit bulk volume if thermal simulation is specified, pressureof one of the fluid phases, elastic rock displacements in the x,y,zdirections and porosity. Some of these variables are eliminated beforeentering the solver. In particular, one of the fluid molar densities iseliminated using the volume balance equation. Usually, a fully implicitsolution of all equations is unnecessary because in many parts of thegrid the fluid throughput in a mass conservation grid cell is not high.In this case an adaptive implicit solution of the equation set is usedin which some of the grid blocks are solved fully implicitly, whileothers are solved using the IMPES (implicit pressure, explicitsaturation) solution scheme. Often the IMPES scheme is used everywherealthough the timestepping of this particular scheme is then limited bythe throughput. The rock force balance equations are always solved fullyimplicitly, even in grid cells which are IMPES.

[0072] The linear solver used to solve the coupled system of non-linearequations is a nested factorization technique²⁵. One such linear solveris disclosed in U.S. Pat. No. 6,230,101 to Wallis, the disclosure ofwhich is incorporated by reference into this specification. Although therock stress equations (1) are nonlinear because of the implicit porosityand highly elliptic and, in addition, are placed next to massconservation equations which display both hyperbolic and parabolic formin the jacobian, the solver has demonstrated good robustness in solvingthe coupled system of equations.

[0073] Although some authors^(8,26) have described fully or partiallycoupled schemes where porosity is calculated directly using rockdilation, ε_(x)+ε_(y)+ε_(z), we have chosen the additional rockconservation equation for two reasons. Firstly, the approach is moregeneral and allows rock to be produced in cases of sand fluidization,wellbore instability, etc. Secondly, as discussed further below, at thebeginning of a new timestep a large plastic displacement of the rock isaccounted for exactly by including it in the implicit rock mass balancefor this timestep.

[0074] Plasticity

[0075] An explicit plastic calculation has been implemented. TheMohr-Coulomb plastic calculation is available. The user can specify thecohesion and angle of internal friction in chosen regions in the grid.At the start of a new timestep, principal stresses are calculated ineach rock force-balance control volume based on the stress field fromthe last timestep. Diameters of the smallest and largest Mohr circlesare calculated and the test is made as to whether failure has occurred.If so, the displacement around which the control volume is centred isaltered to bring the circle back to the failure envelope. Thiscorrective displacement is stored as a plastic displacement. Theseplastic displacements, which can be referred to as (u_(p), v_(p),w_(p)), are part of the total displacements (u, v, w) which in turn area sum of the elastic and plastic contributions.

[0076] Boundary Conditions and Initialisation

[0077] Either stress or displacement boundary conditions can be chosen.Each side of the simulation grid, that is x−, x+, y−, y+ and z−, z+, canhave one of the above boundary conditions prescribed. Usually the bottomof the grid will have a zero displacement condition while the sides andtop will have specified stresses.

[0078] Initialisation of a reservoir with stress boundary conditionsrequires an initial small simulation to equilibrate the simulated field.Before the wells are turned on, boundary stresses are ramped up atintervals until the desired vertical and horizontal stresses areachieved. The pore pressure is maintained at an initial level byincluding an extra production well with a BHP limit set to this level.Compression of the reservoir by the boundary stresses causes the rockmatrix to compress which in turn forces some of the fluid out of thiswell. Fluid-in-place must be checked before starting the simulation.

[0079] User Interface

[0080] The user enters a Young's modulus, Poisson ratio and specificrock density for each grid block. It is also possible to input the fullstiffness matrix, as illustrated above in equation (4), to allow ananisotropic stress tensor. If a plastic calculation is required, theuser enters a cohesion and angle of internal friction in differentregions of the grid.

[0081] Efficiency of Geomechanical Computations

[0082] Because there are several extra equations to be set up andsolved, there is a sizable overhead associated with the geomechanicscalculation. This overhead is offset by several factors. Firstly, thecoupled set of equations is stable and robust. Secondly, implementationof the jacobian and right-hand-side setup includes vectorizationwherever possible. Also, parallel options exist within the simulator.

[0083] On average, for three dimensional IMPES simulations we have notedthat run times are about three times as long as the run time needed forthe same calculation without geomechanics. This is efficient consideringthe IMPES calculation only solves a single equation, whereas fouradditional equations are solved when geomechanics is included. Both theefficiency of setting up the equations and the robustness of the solvercontribute to the overall efficiency. The simulator is capable ofachieving good parallel speedups and this also helps to offset theincreased computational demands.

[0084] Benchmarking Against an Industry Standard Stress Code

[0085] We have chosen the ABAQUS²⁷ stress simulator to test the rockstress calculations in our reservoir simulator. ABAQUS was chosen fortwo reasons. Firstly, it is widely accepted and established. Secondly,it contains a simple, single-phase flow-in-porous-media option.

[0086] Two elastic, one-dimensional problems were designed and run inboth simulators. The test cases differed only in that one featured anaverage sandstone, the other a weak sandstone. They were designed tohave large pore pressure gradients so that the fluid-rock interactionwas pronounced.

[0087] In FIG. 4d, Table 1 presents the simulation parameters used inthe test problems. These are listed in the table as Cases 1 and 2. Awater injection rate of 1000 BBL/D was used. The injector was situatedin one corner, the producer in the far comer. Boundaries were rigid.Dimensions of the grid were chosen as representative of a small patternof injectors and producers. The reservoir simulator was run to steadystate for these cases. In FIG. 4d, under Table 1, in all cases,porosity=0.33, fluid permeability=6 mD, initial pressure=4790 psia,injection/production rate=1000 BBL/D and grid dimensions were 1000ft×1000 ft×20 f.

[0088] In FIG. 2, FIG. 2 shows a comparison of rock displacements aspredicted by the two simulators for Case 1, the average strengthsandstone.

[0089] In FIGS. 3A and 3B, these figures illustrate a comparison for theweak rock. Much larger rock displacements are evident as a result of thevery low Young's modulus. In both cases, the comparison against ABAQUSwas within engineering accuracy.

[0090] Other Test Cases

[0091] In FIG. 4d, Table 1 is illustrated (Table 1 illustratingsimulation parameters used in comparing ABAQUS and Reservoir SimulatorPredictions). We have also run ABAQUS on a further set of 4 problems.Table 1, shown in FIG. 4d, lists these as Cases 3 to 6. These are a pairof two-dimensional and a pair of three-dimensional runs. Forequilibration purposes, the ABAQUS simulations used boundary conditionsthat are not available in our reservoir simulator. We are currentlyworking to remedy this. Initial comparisons in which we haveapproximated the ABAQUS simulation have been encouraging. Bothquantitative and qualitative agreement has been achieved.

[0092] A plastic case was run in ABAQUS to compare against the reservoirsimulator plastic calculation. It is the same as Case 6 in Table 1 ofFIG. 4d. A cohesion and angle of internal friction was input to definethe Mohr-Coulomb failure curve. Differences in boundary conditions inthe two simulators, as discussed above, have prohibited a directcomparison but, again, there is qualitative agreement when comparingplastic displacements. The plastic calculation is robust andtimestepping is not adversely affected.

[0093] In FIGS. 4a 1 and 4 a 2, 4 b 1 and 4 b 2, and 4 c 1 and 4 c 2, wehave simulated some larger fields with the geomechanics option. Anexample is shown in FIGS. 4a 1 and 4 a 2. In this three-dimensional,compositional (10 component) problem, 7 injection and 6 production wellscontribute to the rock stress. Contours of the total normal stress,σ_(x)+σ_(y)+σ_(z) which is also the first stress invariant, are shown inthe mid-xy plane of the simulation. This first stress invariant isimportant because it largely governs compaction processes. Timesteppingwith the geomechanics calculation included was the same as without.There was a factor of 2.8 slowdown in the overall runtime withgeomechanics.

[0094] In summary, the implementation of elastic/plastic stressequations into a commercial, finite difference reservoir simulator hasbeen discussed in this specification. Our goal of demonstratingrobustness together with a comprehensive geomechanical option has beenachieved. Work is continuing to develop this option further. Inparticular:

[0095] 1. The coupled system of equations has demonstrated goodstability and robustness. In some test cases, there was no difference intimestepping when the simulation was run with or without thegeomechanics option.

[0096] 2. Our current capabilities are comprehensive enough to predictboth elastic and plastic rock displacement. The basic design is generalenough to allow future engineering development. This could includewellbore stability and failure, sand production, more accurate faultmodelling.

[0097] 3. We have indicated at various points of this paper our areas ofcurrent work and investigation. In particular these includestress-permeability relationships and modifying our geomechanicalboundary conditions to compare more exactly with ABAQUS. We are alsocurrently working with a large-scale simulation of a very brittlecarbonate field to predict subsidence.

[0098] Nomenclature

[0099] E=Young's modulus

[0100] G=modulus of rigidity, Lamé constant

[0101] H=pressure head, including the gravitational head

[0102] K=fluid phase permeability of rock matrix

[0103] q=single phase fluid flux in a porous medium

[0104] u=rock displacement in the x-direction

[0105] v=rock displacement in the y-direction

[0106] V_(s)=velocity of the solid, as related to Darcy's law

[0107] w=rock displacement in the z-direction

[0108] σ_(x,y,z)=elastic normal rock stress in x,y,z directions

[0109] ε_(x,y,z)=x,y,z elongation strains

[0110] γ_(xy,yz,zx)=shear strains

[0111] λ=Lamé constant

[0112] φ=porosity

[0113] ν=Poisson's ratio

[0114] ρ_(rock)=rock specific density

[0115] τ_(xy,yz,xz)=elastic shear stress

DETAILED DESCRIPTION OF THE INVENTION

[0116] Referring to FIGS. 5 and 6, structured and unstructured grids areillustrated.

[0117] In FIG. 5, an earth formation 15 is illustrated, the formation 15including four (4) horizons 13 which traverse the longitudinal extent ofthe formation 15 in FIG. 5. Recall that a “horizon” 13 is defined to bethe top surface of an earth formation layer, the earth formation layercomprising, for example, sand or shale or limestone, etc.

[0118] A “grid” is located intermediate to the horizon layers 13. Moreparticularly, a “grid” is formed: (1) in between the horizons 13, (2) ontop of the uppermost horizon 13, and (3) below the lowermost horizon 13in the Earth formation 15. When gridding the formation 15, the act orfunction of “gridding” the formation 15 includes the step of dividingthe formation 15 into a multitude of individual “cells” which, whenconnected together, comprise the “grid”.

[0119] In FIG. 6, for example, the Earth formation 15 includes anuppermost horizon 13 a and a lowermost horizon 13 b which is separatedfrom the uppermost horizon 13 a by an intermediate earth formation layer15 a. The intermediate earth formation layer 15 a includes, for example,a sand layer or a shale layer or a limestone layer, etc. A particular“gridding software” will “grid” the earth formation layer 15 a; that is,the formation layer 15 a will be divided up, by the “gridding software”into a multitude of cells 15 a 1.

[0120] In the prior art, a “gridding software” product known as “Grid”was marketed by GeoQuest, a division of Schlumberger TechnologyCorporation, Abingdon, the United Kingdom (UK). The “Grid” softwarewould divide the formation layers 15 a into a multitude of cells.However, each of the multitude of cells were approximately “rectangular”in cross sectional shape.

[0121] In addition, another “gridding software”, known as “Petragrid” isdisclosed in U.S. Pat. No. 6,018,497, the disclosure of which isincorporated by reference into this specification. The “Petragrid”gridding software will ‘grid’ the formation withtriangularly/tetrahedrally shaped grid cells, known as “unstructured ridcells”.

[0122] In addition, another “gridding software”, known as “Flogrid”, isdisclosed in U.S. Pat. No. 6,106,561, the disclosure of which isincorporated by reference into this specification. The “Flogrid”gridding software includes the “Petragrid” gridding software; however,in addition, the “Flogrid” gridding software will ‘grid’ the formationwith rectangularly shaped grid cells, known as “structured grid cells”.

[0123] In FIG. 6, the cells 15 a 1 are shown to be approximately“rectangular” in cross sectional shape. These grid cells are “structuredgrid cells”.

[0124] In FIG. 5, however, a multitude of grid cells 15 a 1 have beenformed in the earth formation 15 intermediate the horizons 13, and eachcell 15 a 1 has a cross sectional shape that, in addition to beingapproximately “rectangular” in cross section, is either approximately“polygonal” or “tetrahedral” in cross section. FIG. 5 clearly shows amultitude of cells 15 a 1 where each cell 15 a 1 has a cross sectionalshape which is either approximately “polygonal” or “tetrahedral” incross sectional shape (i.e., an “unstructured grid cell”) in addition tobeing approximately “rectangular” in cross sectional shape (i.e., a“structured grid cell”).

[0125] Referring to FIGS. 7, 8, and 9, other examples of “structured”and “unstructured” grid cells are illustrated. FIG. 7 illustratesexamples of “un-structured” grid cells having a triangular/tetrahedralcross sectional shape. FIG. 8 illustrates examples of “structured” gridcells having an approximate “rectangular” cross sectional shape. FIG. 9illustrates further examples of “unstructured” grid cells having atriangular/tetrahedral cross sectional shape.

[0126] Referring to FIGS. 10 and 11a, a “staggered grid” will now bedefined with reference to FIGS. 10 and 11a.

[0127] In FIG. 10, a plurality of “structured” grid cells areillustrated (i.e., cells having an approximately rectangularly shapedcross sectional shape). However, in FIG. 10, a pair of “staggered gridcells” are also illustrated. Assuming that a ‘first structured orunstructured grid cell’ is disposed directly adjacent a ‘secondneighboring structured or unstructured grid cell’, a “staggered gridcell” can be defined as one which consists of a ‘first half’ and a‘second half’, where the ‘first half’ of the ‘staggered grid cell’comprises one-half of the ‘first structured or unstructured grid cell’and the ‘second half’ of the ‘staggered grid cell’ comprises one-half ofthe ‘second neighboring structured or unstructured grid cell’.

[0128] For example, in FIG. 10, the staggered grid cell 20 includes theone-half 20 a of a first structured grid cell and the one-half 20 b ofthe neighboring, second structured grid cell. Similarly, the staggeredgrid cell 22 includes the one-half 22 a of a first structured grid celland the one-half 22 b of a neighboring, second structured grid cell.

[0129] In FIG. 11a, a pair of “unstructured” grid cells are illustrated(i.e., cells each having an approximately triangularly or tetrahedrallyshaped cross sectional shape). In FIG. 11a, a staggered grid cell 24consists of a first half 24 a from one unstructured grid cell 26 and asecond half 24 b from a second unstructured grid cell 28.

[0130] Referring to FIGS. 11b, 12, 13, and 14, a computer workstation 44or other computer system 44 is illustrated.

[0131] In FIG. 11b, the computer workstation 44 in FIG. 11b includes asystem bus, a workstation processor 44 d connected to the bus, aworkstation memory 44 a connected to the bus, and a recorder or displayor 3D viewer 44 e connected to the bus. A CD-Rom 46 is inserted into theworkstation and the following software is loaded from the CD-Rom 46 andinto the workstation memory 44 a: (1) a “Flogrid” software 46 a, and (2)an Eclipse simulator software 46 b. The “Flogrid” software is describedin U.S. Pat. No. 6,106,561 to Farmer, the disclosure of which hasalready been incorporated by reference into this specification. Inputdata is provided to and input to the workstation 44: (1) a well logoutput record 42, and (2) a reduced seismic data output record 24 a.

[0132] In FIG. 12, the workstation memory 44 a of FIG. 11b stores theFlogrid software 46 a and the Eclipse simulator software 46 b. TheFlogrid software 46 a includes a reservoir data store, a reservoirframework, a structured gridder, a “Petragrid” Un-structured gridder,and an Upscaler. The outputs 47 from the Upscaler and from the PetragridUnstructured gridder are provided as inputs to the Eclipse simulatorsoftware 46 b. The “Petragrid” Un-structured gridder is disclosed inU.S. Pat. No. 6,018,497 to Gunasekera, the disclosure of which hasalready been incorporated by reference into this specification. TheEclipse simulator software 46 b, responsive to the outputs 47, generatesa set of ‘simulation results’ 48 which are provided to the 3D viewer 44e. The Eclipse simulator software 46 b includes a “Linear Solver” whichis discussed in U.S. Pat. No. 6,230,101 to Wallis, the disclosure ofwhich is incorporated by reference into this specification.

[0133] In FIG. 13, to summarize FIGS. 11b and 12, output data 47 isgenerated from the Upscaler and the Petragrid un-structured gridder inthe Flogrid software 46 a, and that output data 47 is provided to theEclipse simulator software 46 b, along with the outputs of otherprograms 49. Responsive thereto, the Eclipse simulator software 46 bgenerates a set of simulation results 48 which are, in turn, provided tothe 3D viewer 44 e.

[0134] In FIG. 14, one example of the simulation results 48, which aredisplayed on the 3D viewer 44 e of FIG. 13, is illustrated.

[0135] Referring to FIG. 15, as previously discussed, the workstationmemory 44 a stores the Eclipse simulator software 46 b. However, inaccordance with one aspect of the present invention, the Eclipsesimulator software 46 b includes a Parameter Determination software 50.The “parameters” are determined for each “staggered grid cell”, such aseach of the “staggered grid cells” shown in FIGS. 10 and 11a of thedrawings. The “parameters”, which are determined by the ParameterDetermination software 50 for each ‘staggered grid cell’, include thefollowing “parameters”, which have been previously discussed in the“Description of the Invention” portion of this specification:

[0136] u=rock displacement in the x-direction

[0137] v=rock displacement in the y-direction

[0138] w=rock displacement in the z-direction

[0139] Therefore, the parameter (u, v, w) represent the movement of therock. When the aforementioned parameters (u, v, w) have been determined,the following “additional parameters” are determined for each ‘staggeredgrid cell’ from and in response to the parameters (u, v, w), and these“additional parameters” have also been discussed in the “Description ofthe Invention” portion of this specification:

[0140] ε_(x,y,z)=x,y,z elongation strains

[0141] γ_(xy,yz,zx)=shear strains

[0142] σ_(x,y,z)=elastic normal rock stress in x,y,z directions

[0143] τ_(xy,yz,xz)=elastic shear stress

[0144] φ=porosity of rock (variable)

[0145] u_(p)=plastic rock displacement in the x-direction

[0146] v_(p)=plastic rock displacement in the y-direction

[0147] w_(p)=plastic rock displacement in the z-direction

[0148] Recall equation (3) set forth in the “Description of theInvention” portion of this specification, as follows: $\begin{matrix}{\begin{matrix}{ɛ_{x} = \frac{\partial u}{\partial x}} \\{ɛ_{y} = \frac{\partial v}{\partial y}} \\{ɛ_{z} = \frac{\partial w}{\partial z}} \\{\gamma_{xy} = {\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}}} \\{\gamma_{yz} = {\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}}} \\{\gamma_{zx} = {\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}}}\end{matrix}.} & (3)\end{matrix}$

[0149] We start by estimating values for (u, v, w), which is hereinafterreferred to as the ‘estimating procedure’. Therefore, when (u, v, w) isestimated by the Parameter Determination software 50 of the presentinvention, the values of “ε_(x,y,z)” (the ‘x,y,z elongation strains’)and “γ_(xy,yz,zx)” (the ‘shear strains’) are determined from the aboveequation (3). When “ε_(x,y,z)” and “γ_(xy,yz,zx)” are determined fromequation (3), the values of “σ_(x,y,z)” (the ‘elastic normal rock stressin x,y,z directions’) and “τ_(xy,yz,xz)” (the ‘elastic shear stress’)are further determined from the above described equation (2), asfollows: $\begin{matrix}{\begin{matrix}{\sigma_{x} = {{2G\quad ɛ_{x}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\sigma_{y} = {{2G\quad ɛ_{y}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\sigma_{z} = {{2G\quad ɛ_{z}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\tau_{xy} = {G\quad \gamma_{xy}}} \\{\tau_{yz} = {G\quad \gamma_{yz}}} \\{\tau_{zx} = G_{\gamma_{zx}}}\end{matrix}.} & (2)\end{matrix}$

[0150] Now that “σ_(x,y,z)” and “τ_(xy,yz,xz)” are known, the previouslydescribed set of ‘rock momentum balance’ differential equations [thatis, equation (1)] are solved using the aforementioned known values of“σ_(x,y,z)” and “τ_(xy,yz,xz)”. After the ‘rock momentum balance’differential equations are solved, if the resultant ‘residuals’ aredetermined to be approximately equal to zero (0), the previouslyestimated values of (u, v, w) are thereby determined to represent‘accurate rock displacement parameters’ for the reservoir. Subsidence isdetermined from the ‘accurate rock displacement parameters’ by solving afinal set of failure criterion equations, which comprise the residualsand any derivatives, to determine a set of rock plastic displacements(u_(p)=plastic rock displacement in the x-direction, v_(p)=plastic rockdisplacement in the y-direction, and w_(p)=plastic rock displacement inthe z-direction) forming a part of the rock displacement parameters (u,v, w).

[0151] The Parameter Determination software 50 will also determine “φ”,the porosity of rock, which is a variable set forth in equation (5) setforth above.

[0152] Referring to FIG. 16, when the parameters (u, v, w, ε_(x,y,z),γ_(xy,yz,zx), σ_(x,y,z), τ_(xy,yz,xz), φ, u_(p), v_(p), w_(p)) have beendetermined, as discussed above, for each staggered grid block or cell,the Recorder or Display or 3D Viewer 44 e of FIG. 11b will illustrate a“Table” similar to the “Table” shown in FIG. 16.

[0153] In the Table of FIG. 16, for staggered grid block 1 (“SGB1”), afirst set of the parameters “(u, v, w, ε_(x,y,z), γ_(xy,yz,zx),σ_(x,y,z), τ_(xy,yz,xz), φ, u_(p), v_(p), w_(p))” will be generated anddisplayed on the Recorder or Display device 44 e.

[0154] Similarly, in the Table of FIG. 16, for staggered grid block 2(“SGB2”), a second set of the parameters “(u, v, w, ε_(x,y,z),γ_(xy,yz,zx), σ_(x,y,z), τ_(xy,yz,xz), φ, u_(p), v_(p), w_(p))” will bedisplayed on the Recorder or Display device 44 e.

[0155] Similarly, in the Table of FIG. 16, for staggered grid block “n”(“SGBn”), an n-th set of the parameters “(u, v, w, ε_(x,y,z),γ_(xy,yz,zx), σ_(x,y,z), τ_(xy,yz,xz), φ, u_(p), v_(p), w_(p))” will begenerated and displayed on the Recorder or Display device 44 e.

[0156] Referring to FIG. 17, when the “Table” of FIG. 16 is generatedand displayed on the display device 44 e, it is now possible todetermine ‘Subsidence’ in an oilfield reservoir. The term ‘subsidence’refers to a failing or ‘giving away’ of the seabed floor. In FIG. 17, adrilling rig 52 is situated at the surface of the sea 54. The drillingrig 52 withdraws underground deposits of hydrocarbon (e.g., oil) andother fluids (e.g. water) from an undersea reservoir 56 which is locatedbelow the seabed floor 58. Over a period of, for example, ten years, acertain amount of subsidence 60 occurs due to the withdrawal of theunderground deposits of hydrocarbon and water from the underseareservoir 56, the subsidence 60 producing a lowering of the seabed floor58 a certain amount which corresponds to the subsidence 60.

[0157] In FIG. 17, the Parameter Determination Software 50 of FIG. 15determines the aforementioned subsidence 60 by first determining (viathe ‘estimating procedure’ discussed above) the ‘accurate rockdisplacement parameters (u, v, w)’ for the reservoir. These ‘accuraterock displacement parameters (u, v, w)’ represent the ‘movement of therock’, or, in our example, the displacements of the seabed floor 58 inthe (x, y, z) directions over the ten year example period. Thesedisplacements of the seabed floor 58 occur as a result of the subsidence60 (that is, the ‘failing’ or the ‘giving away’ of the ground or seabedfloor) illustrated in FIG. 17. Therefore, the parameters (u, v, w)represent or characterize ‘rock movement in the x, y, and z directions’which, in turn, represent or characterize the ‘subsidence 60’. Note thatthe ‘accurate rock displacement parameters (u, v, w)’ may include: (1) adisplacement that is ‘elastic’, and (2) a displacement that is ‘plastic’which is denoted (u_(p), v_(p), w_(p)). These plastic displacements(u_(p), v_(p), w_(p)) were referred to above under ‘Description of theInvention’ in the section entitled ‘Plasticity’. ‘Subsidence’ is theresult of a ‘failure’ of the rock, for example by crushing, cracking orsome other failure mechanism. When the rock has failed, some of itsdisplacement is not recoverable when the original conditions areimposed, and it is the presence of this unrecoverable displacement(u_(p), v_(p), w_(p)) that characterizes the ‘subsidence’ 60. Also notethat the undersea formation in FIG. 17 is gridded with structured andunstructured grids, and, as a result, the undersea formation in FIG. 17is gridded with “staggered grids”, as graphically illustrated by thestaggered grid 62 in FIG. 17. The flowchart of FIG. 18 will describe themethod, practiced by the Parameter Determination Software 50 of FIG. 15,for determining the ‘accurate rock displacement parameters (u, v, w)’which represent these displacements in the (x, y, z) directions due tothe subsidence 60.

[0158] Referring to FIG. 18, a flowchart is illustrated which describesthe method practiced by the Parameter Determination Software 50 fordetermining the parameters (u, v, w) representing the rock displacementsin the (x, y, z) directions (i.e., the ‘movement of the rock’) due tothe subsidence illustrated in FIG. 17.

[0159] In FIG. 18, the Parameter Determination Software 50 willdetermine the ‘rock movement’ parameters (u, v, w) in accordance withthe following steps. With reference to FIG. 18, each step of theParameter Determination software 50 will be discussed in detail below asfollows:

[0160] 1. Input User Defined Rock Mechanical Properties—Block 50 a

[0161] The user inputs a lot of data, called the “input file”,comprising keywords that define the particular simulation that isdesired. Part of these keywords represent the mechanical properties ofthe rock. These “mechanical properties of the rock” are the kind ofproperties that are needed to drive the equations set forth above knownas the “rock momentum balance equations” which are identified above as“equation (1)”.

[0162] 2. Build Staggered Grid Volumes, Areas into GeomechanicalArrays—Block 50 b

[0163] When the “input file” comprising the “mechanical properties ofthe rock” is entered, the staggered grid can be built. In this case, weneed the staggered grid ‘volumes’ (how much volume space the staggeredgrid would enclose) and the ‘areas’ (how much surface area the staggeredgrids would include). These ‘volumes’ and ‘areas’ are built into special‘arrays’ used for geomechanics calculations.

[0164] 3. Build Special NNC's for Staggered Grid in NNC Structure—Block50 c

[0165] When the aforementioned ‘arrays’, used for the geomechanicscalculations, are built, the next step is to build special NNC's. Theacronym NNC refers to “Non-Neighbor Connections”. Because the abovedefined ‘equation (1)’ [the “rock momentum balance equations”] are moredifficult to solve than the standard reservoir simulation equations, itis necessary to create extra connections on the different grid blocksduring the simulation in order to solve ‘equation (1)’ properly.Therefore, during this step, we build up these “Non-NeighborConnections” for the staggered grid in a special “Non-NeighborConnections Structure” or “NNC structure” that has been created insidethe simulator. The “NNC structure” is a collection of arrays and specialareas that have been set aside inside the simulator.

[0166] 4. Build Boundary Structure Based on Staggered Gridding—Block 50d

[0167] When the “NNC structure” has been created inside the simulator,it is necessary to build up the boundaries representing the edges of thesimulation or grid. The edges of the grid need this ‘boundarystructure’, which is, in turn, based on the staggered grid. The whole‘boundary structure’, which is a collection or arrays representing the‘boundary conditions’, has to be set up in order to rigorously solve theabove mentioned equations, namely, the “rock momentum balance equations”defined above as ‘equation (1)’.

[0168] 5. Assign Input Mechanical Parameters to Staggered Grid ArrayStructure—Block 50 e

[0169] When the ‘boundary structure’ has been build based on thestaggered gridding, it is necessary to assign input mechanicalparameters to the staggered grid array structure. After this isaccomplished, the time stepping can now commence.

[0170] 6. Start Time Stepping—Block 50 f

[0171] When the input mechanical parameters are assigned to thestaggered grid array structure, the simulator will now start ‘timestepping’. The simulator 46 b of FIG. 15 will time step forward in timeand project into the future to determine what will happen.

[0172] 7. Calculate and Store Equation Residuals for Rock and FluidForce Balance in Staggered Grid—Block 50 g

[0173] The first thing that is accomplished after the ‘time stepping’ iscommenced is to calculate and store the equation residuals for the ‘rockand fluid force balance’. The ‘rock and fluid force balance’ refers to‘equation (1)’, that is, the “rock momentum balance equations”, whichare set forth above and are duplicated below as follows: $\begin{matrix}{\begin{matrix}{{\frac{\partial\sigma_{x}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{x} + P_{x}} = 0} \\{{\frac{\partial\sigma_{y}}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{zy}}{\partial z} + F_{y} + P_{y}} = 0} \\{{\frac{\partial\sigma_{x}}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + F_{z} + P_{z}} = 0}\end{matrix}.} & (1)\end{matrix}$

[0174] The simulator will actually calculate the “rock momentum balanceequations” representing ‘equation (1)’, and the residuals are exactly asset forth in ‘equation (1)’. Referring to ‘equation (1)’ set forth above(that is, the “rock momentum balance equations”), the term “residual”can be defined as follows: “How different from ‘0’ is the left-hand sideof the above referenced ‘equation (1)’?”. The objective is to try tomake the left-hand side of ‘equation (1)’ equal to ‘0’, which is theright-hand side of the ‘equation (1)’. The simulator 46 b of FIG. 15will therefore “calculate the ‘rock momentum balance equations’representing ‘equation (1)’, that is, the simulator 46 b will adjust theparameters (u, v, w) and porosity and other resultant variables as setforth above in equations (2) and (3), until the left-hand side of theabove ‘equation (1)’ is equal to the right-hand side of ‘equation (1)’,where the right-hand side of ‘equation (1)’ is equal to ‘0’. Therefore,we know that we have solved the above referenced “rock momentum balanceequations” of ‘equation (1)’ when the left-hand side of ‘equation (1)’is equal to ‘0’. We find out how far away from ‘0’ we are by calculatingthe ‘residuals’.

[0175] 8. Calculate and Store Equation Residual Derivatives for Rock andFluid Force Balance—Block 50 h

[0176] At this point, we calculate and store the equation residualderivatives (which are required for the Newton gradient search), whichderivatives will drive the “rock momentum balance equations” residuals[identified as ‘equation (1)’] to ‘0’. At this point, we have: (1)‘residuals’, and we have (2) ‘derivatives of these residuals’.

[0177] 9. Solve This Linear System—Block 50 i

[0178] Now that we know:

[0179] (1) the ‘Residuals’, and

[0180] (2) the ‘Derivatives of these Residuals’,

[0181] we can now solve this whole ‘system of linear equations’representing the ‘rock momentum balance equations’ of ‘equation (1)’which are set forth again below, as follows: $\begin{matrix}{\begin{matrix}{{\frac{\partial\sigma_{x}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{x} + P_{x}} = 0} \\{{\frac{\partial\sigma_{y}}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{zy}}{\partial z} + F_{y} + P_{y}} = 0} \\{{\frac{\partial\sigma_{x}}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + F_{z} + P_{z}} = 0}\end{matrix}.} & (1)\end{matrix}$

[0182] The above referenced ‘rock momentum balance equations’ of‘equation (1)’ are solved together and simultaneously with the standardreservoir simulation equations.

[0183] When the above referenced ‘rock momentum balance equations’ of‘equation (1)’ are solved, the “u”, “v”, and “w” displacementparameters, representing movement of the rock, are determined, where:

[0184] u=rock displacement in the x-direction

[0185] v=rock displacement in the y-direction

[0186] w=rock displacement in the z-direction

[0187] When the displacement parameters (u, v, w) are determined, theabove referenced ‘equation (3)’ is used to determine “ε_(x,y,z)” (the‘x,y,z elongation strains’) and “γ_(xy,yz,zx)” (the ‘shear strains’),since “ε_(x,y,z)” and “γ_(xy,yz,zx)” are function of “u”, “v”, and “w”,as follows: $\begin{matrix}\begin{matrix}{ɛ_{x} = \frac{\partial u}{\partial x}} \\{ɛ_{y} = \frac{\partial v}{\partial y}} \\{ɛ_{z} = \frac{\partial w}{\partial z}} \\{\gamma_{xy} = {\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}}} \\{\gamma_{yz} = {\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}}} \\{\gamma_{zx} = {\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}}}\end{matrix} & (3)\end{matrix}$

[0188] When “ε_(x,y,z)” and “γ_(xy,yz,zx)” are determined from ‘equation(3)’, the above referenced ‘equation (2)’ is used to determine“σ_(x,y,z)” (the ‘elastic normal rock stress in x,y,z directions’) and“τ_(xy,yz,xz)” (the ‘elastic shear stress’) since “σ_(x,y,z)” and“τ_(xy,yz,xz)” are a function of “ε_(x,y,z)” and “γ_(xy,yz,zx)” in‘equation (2)’ as follows: $\begin{matrix}\begin{matrix}{\sigma_{x} = {{2G\quad ɛ_{x}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\sigma_{y} = {{2G\quad ɛ_{y}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\sigma_{z} = {{2G\quad ɛ_{z}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\{\tau_{xy} = {G\quad \gamma_{xy}}} \\{\tau_{yz} = {G\quad \gamma_{yz}}} \\{\tau_{zx} = G_{\gamma_{zx}}}\end{matrix} & (2)\end{matrix}$

[0189] When the “σ_(x,y,z)” and “τ_(xy,yz,xz)” are determined from‘equation (2)’, the above referenced ‘equation (1)’ is solved, since‘equation (1)’ is a function of “σ_(x,y,z)” and “τ_(xy,yz,xz)”, asfollows: $\begin{matrix}\begin{matrix}{{\frac{\partial\sigma_{x}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{x} + P_{x}} = 0} \\{{\frac{\partial\sigma_{y}}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{zy}}{\partial z} + F_{y} + P_{y}} = 0} \\{{\frac{\partial\sigma_{x}}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + F_{z} + P_{z}} = 0}\end{matrix} & (1)\end{matrix}$

[0190] The result of these calculations will produce the Table of FIG.16. A reservoir engineer would be interested in knowing each of thequantities shown in the Table of FIG. 16.

[0191] 10. Evaluate Failure Criteria and Solve for PlasticDisplacements—Block 50 j

[0192] A final step is required if the reservoir engineer has specifiedregions in the step 50 a in FIG. 18 where the rock may fail. In thiscase, the final “σ_(x,y,z)” and “τ_(xy,yz,xz)” are used to evaluate thechosen failure criterion which may be Mohr-Coulomb as referred to abovein the section ‘Description of the Invention’ subsection Plasticity. Ifthe failure criterion is exceeded, then a set of equations is set up tosolve for the plastic displacements, denoted (u_(p), v_(p), w_(p))above, which will force the failure criterion residual to be exactlyzero. This system is set up and solved exactly as described above withthe exception that the residuals are now the failure criterion residualswhich are driven to zero instead of the equation (1) residuals. Thisstep is 50 j in FIG. 18.

[0193] In FIG. 18, when the ‘system of linear equations’ set forth aboveas ‘equation (1)’ has been solved, referring to the ‘feed-back loop’ 64in FIG. 18, the simulator 46 b of FIG. 15 continues its ‘time stepping’by repeating the implementation of blocks 50 g, 50 h, 50 i and 50 j inFIG. 18.

[0194] When the flow chart of FIG. 18 is executed, the ‘result’ is: asolution to a problem, which a reservoir engineer would want to solve,in which you have certain boundary conditions applied to the rock,boundary conditions applied to ‘equation (1)’, and boundary conditionsapplied to the reservoir simulation equations. For example, in the NorthSea, problems with ‘subsidence’ exist. The term ‘subsidence’ is asituation where the ground actually ‘gives-way’ or ‘fails’. The abovecalculation would be applied to this ‘subsidence’ problem. The abovetype of calculation, which represents a solution to the above referenced‘Elastic Stress Equations’ using ‘staggered gridding’, would be used tohandle the ‘more difficult’ fields, such as the fields having the‘subsidence’ problem of FIG. 17.

[0195] In summary, rock will initially behave elastically because it isin equilibrium and there are no external forces. Elastic behavior meansthat the rock behaves like a spring; that is, you can push it (i.e.,apply force) and it will return exactly to its original position whenyou stop pushing. Subsidence happens when the rock is pushed beyond itsown strength. For example, in a particular reservoir where certain wellsare withdrawing fluids, the rock may no longer have the support of thewater pressure, and the weight of the overburden will become high enoughto cause the rock to fail. If the rock fails, e.g. by crushing orcracking, then, it no longer moves elastically but rather movesplastically. This plastic movement is what is referred to as“subsidence”. We know this will occur because we solve ‘equation (1)’,which are the rock momentum balance equations (otherwise known as rockforce balance equations). If these equations tell us that the forcesbecome too large in one direction relative to the forces in anotherdirection at some particular point in time, then, we know that failurewill occur. At this point, we further solve a final set of equationswhich will calculate the plastic displacements. “Derivatives” are neededbecause these equations are nonlinear.

[0196] References

[0197] The following references (1 through 26) are incorporated byreference into the specification of this application:

[0198] 1. Hansen, K. S., Prats, M., Chan, C. K., “Finite-elementmodeling of depletion-induced reservoir compaction and surfacesubsidence in the South Belridge oil field, California”, SPE 26074,presented at the Western Regional Meeting, Anchorage, Ak., U.S.A., May26-28, 1993.

[0199] 2. Berumen, S., Cipolla, C., Finkbeiner, T., Wolhart, S.,Rodriguez, F., “Integrated reservoir geomechanics techniques in theBurgos Basin, Mexico: An improved gas reservoirs management”, SPE 59418presented at the 2000 SPE Asia Pacific Conference on IntegratedModelling for Asset Management, Yokohama, Japan, Apr. 25-26, 2000.

[0200] 3. Kojic, M., Cheatham, J. B., “Theory of plasticity of porousmedia with fluid flow”, SPE 4394, printed in Transactions, volume 257,June 1974.

[0201] 4. Kojic, M., Cheatham, J. B., “Analysis of the influence offluid flow on the plasticity of porous rock under an axially symmetricpunch”, SPE 4243, printed in Transactions, volume 257, June 1974.

[0202] 5. Corapcioglu, M. Y., Bear, J., “Land Subsidence”, presented atthe NATO Advanced Study Institute on Mechanics of Fluids in PorousMedia, Newark Del., U.S.A., Jul. 18-27, 1982.

[0203] 6. Demirdzic, I., Martinovic, D., “Finite volume method forthermo-elasto-plastic stress analysis”, Computer Method in AppliedMechanics and Engineering 109, 331-349, 1993.

[0204] 7. Demirdzic, I., Muzaferija, S., “Finite volume method forstress analysis on complex domains”, International Journal for NumericalMethods in Engineering, vol. 37, 3751-3766, 1994.

[0205] 8. Settari, A., Walters, D. A., “Advances in coupledgeomechanical and reservoir modeling with applications to reservoircompaction”, SPE 51927 presented at the 1999 SPE Reservoir SimulationSymposium, Houston Tex., Feb. 14-17, 1999.

[0206] 9. Heffer, K. J., Koutsabeloulis, N. C., Wong, S. K., “Coupledgeomechanical, thermal and fluid flow modelling as an aid to improvingwaterflood sweep efficiency”, Eurock '94, Balkema, Rotterdam, 1994.

[0207] 10. Gutierrez, M., and Lewis, R. W., “The role of geomechanics inreservoir simulation”, Proceedings Eurock '98, Vol. 2, 439-448, 1998.

[0208] 11. Geertsma, J., “Land subsidence above compacting oil and gasreservoirs”, JPT (1973), 734-744.

[0209] 12. Holt, R. E., “Permeability reduction induced by anonhydrostatic stress field”, SPE Formation Evauation, 5, 444-448.

[0210] 13. Rhett, W. and Teufel, L. W., “Effect of reservoir stress pathon compressibility and permeability of sandstones”, presented at SPE67^(th) Annual Technical conference and Exhibition, SPE 24756,Washington, D.C., October, 1992.

[0211] 14. Ferfera, F. R., Sarda, J. P., Bouteca, M., and Vincke, O.,“Experimental study of monophasic permeability changes under variousstress paths”, Proceedings of 36^(th) U. S. Rock Mechanics Symposium,Kim (ed.), Elsevier, Amsterdam.

[0212] 15. Teufel, L. W., Rhett, W. and Farrell, H. E., “Effect ofreservoir depletion and pore pressure drawdown on in situ stress anddeformation in the Ekofisk field, North Sea”, Rock Mechanics as aMultidisciplinary Science, Roegiers (ed.), Balkema, Rotterdam.

[0213] 16. Teufel, L. W., and Rhett, W., “Geomechanical evidence forshear failure of chalk during production of the Ekofisk field”, SPE22755 presented at the 1991 SPE 66^(th) Annual Technical Conference andExhibition, Dallas, October 6-9.

[0214] 17. Chou, C. C., Pagano, N. J., “Elasticity, Tensor, Dyadic, andEngineering Approaches”, Dover Publications, Inc., New York, 1967.

[0215] 18. Budynas, R. G., “Advanced Strength and Applied StressAnalysis: Second Edition”, McGraw-Hill Book Co., 1999.

[0216] 19. Naccache, P. F. N., “A fully implicit thermal simulator”, SPE37985 presented at the SPE Reservoir Simulation Symposium, September1997.

[0217] 20. Aziz, K., Settari, A., “Petroleum Reservoir Simulation”,Applied Science Publ., New York City (1979).

[0218] 21. Terzaghi, K., “Erdbaumechanik auf bodenphysikalischerGrundlage”, Franz Deuticke, Vienna, 1925.

[0219] 22. Davies, J. P., Davies, D. K., “Stress-dependent permeability:characterization and modeling”, SPE 56813 presented at the 1999 SPEAnnual Technical Conference and Exhibition held in Houston, Tex., Oct.3-6, 1999.

[0220] 23. Yale, D. P., Crawford, B., “Plasticity and permeability incarbonates: dependence on stress path and porosity”, SPE/ISRM 47582,presented at the SPE/ISRM Eurock '98 Meeting, Trondheim, Norway, Jul.8-10, 1998.

[0221] 24. Appleyard, J. R., Cheshire, I. M., “Nested factorization”,SPE 12264 presented at the SPE Reservoir Simulation Symposium, SanFrancisco, Nov. 15-18, 1983.

[0222] 25. Settari, A., Mourits, F. M., “A coupled reservoir andgeomechanical simulation system”, SPE 50939, SPE Journal, 219-226,September, 1998.

[0223] 26. ABAQUS Theory Manual, Hibbitt, Karlsson & Sorensen, Inc.,1995.

[0224] The invention being thus described, it will be obvious that thesame may be varied in many ways. Such variations are not to be regardedas a departure from the spirit and scope of the invention, and all suchmodifications as would be obvious to one skilled in the art are intendedto be included within the scope of the following claims.

1. A method of determining subsidence in a reservoir, said subsidence resulting from a continual withdrawal of hydrocarbon deposits and other fluids from said reservoir, comprising the steps of: (a) estimating rock displacement parameters u, v, and w representing rock movement in the x, y, and z directions; (b) determining ε_(x,y,z) (‘x,y,z elongation strains’) and γ_(xy,yz,zx) (‘shear strains’) from the rock displacement parameters u, v, and w; (c) determining σ_(x,y,z) (‘elastic normal rock stress in x,y,z directions’) and τ_(xy,yz,xz) (‘elastic shear stress’) from the ε_(x,y,z) and γ_(xy,yz,zx); (d) solving a set of rock momentum balance differential equations from the σ_(x,y,z) and the τ_(xy,yz,xz) and determining if any residuals exist, the rock displacement parameters u, v, and w representing accurate rock displacement parameters for said reservoir when the residuals are substantially equal to zero; and (e) determining said subsidence in said reservoir from said accurate rock displacement parameters determined from step (d), by solving a final set of failure criterion equations, which comprise the residuals and any derivatives, to determine a set of rock plastic displacements, said rock plastic displacements forming a part of said rock displacement parameters u, v, and w.
 2. The method of claim 1, wherein the determining step (b) comprises the step of: determining said ε_(x,y,z) (the ‘x,y,z elongation strains’) and said γ_(xy,yz,zx) (the ‘shear strains’) from said rock displacement parameters u, v, and w by using the following equations: $\begin{matrix} {ɛ_{x} = \frac{\partial u}{\partial x}} \\ {ɛ_{y} = \frac{\partial v}{\partial y}} \\ {ɛ_{z} = \frac{\partial w}{\partial z}} \\ {\gamma_{xy} = {\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}}} \\ {\gamma_{yz} = {\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}}} \\ {\gamma_{zx} = {\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}}} \end{matrix}$


3. The method of claim 2, wherein the determining step (c) comprises the step of: determining said σ_(x,y,z) (the ‘elastic normal rock stress in x,y,z directions’) and said τ_(xy,yz,xz) (the ‘elastic shear stress’) from said ε_(x,y,z) and γ_(xy,yz,zx) by using the following equations: $\begin{matrix} {\sigma_{x} = {{2G\quad ɛ_{x}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\ {\sigma_{y} = {{2G\quad ɛ_{y}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\ {\sigma_{z} = {{2G\quad ɛ_{z}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\ {\tau_{xy} = {G\quad \gamma_{xy}}} \\ {\tau_{yz} = {G\quad \gamma_{yz}}} \\ {\tau_{zx} = G_{\gamma_{zx}}} \end{matrix}$


4. The method of claim 3, wherein the solving step (d) comprises the steps of: solving said set of rock momentum balance differential equations from said σ_(x,y,z) and said τ_(xy,yz,xz), wherein said rock momentum balance differential equations have a left-hand side on a left side of an equal sign and a right-hand side on a right side of said equal sign, said rock momentum balance differential equations including, $\begin{matrix} {{\frac{\partial\sigma_{x}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{x} + P_{x}} = 0} \\ {{\frac{\partial\sigma_{y}}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{zy}}{\partial z} + F_{y} + P_{y}} = 0} \\ {{\frac{\partial\sigma_{z}}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + F_{z} + P_{z}} = 0} \end{matrix}$

said rock displacement parameters u, v, and w estimated in step (a) representing said accurate rock displacement parameters for said reservoir when a left-hand side of said rock momentum balance differential equations substantially equal a right hand side of said rock momentum balance differential equations,
 5. The method of claim 4, wherein the solving step (e) comprises the steps of using the said σ_(x,y,z) and said τ_(xy,yz,xz) to determine whether a failure criterion has been exceeded in one of a set of staggered grids.
 6. The method of claim 5, wherein, when said failure criterion has been exceeded in said one of said set of staggered grids, the solving step (e) further comprises the step of: setting up a final set of equations, each equation being a failure criterion approximately equal to zero, and determining a set of parameters up, vp, wp which are plastic displacements resulting from a failure of the rock, said plastic displacements forming a part of said subsidence and being determined from said accurate rock displacement parameters for said reservoir.
 7. A program storage device readable by a machine, tangibly embodying a program of instructions executable by the machine to perform method steps for determining subsidence in a reservoir, said subsidence resulting from a continual withdrawal of hydrocarbon deposits and other fluids from said reservoir, said method steps comprising: (a) estimating rock displacement parameters u, v, and w representing rock movement in the x, y, and z directions; (b) determining ε_(x,y,z) (‘x,y,z elongation strains’) and γ_(xy,yz,zx) (‘shear strains’) from the rock displacement parameters u, v, and w; (c) determining σ_(x,y,z) (‘elastic normal rock stress in x,y,z directions’) and τ_(xy,yz,xz) (‘elastic shear stress’) from the ε_(x,y,z) and γ_(xy,yz,zx); (d) solving a set of rock momentum balance differential equations from the σ_(x,y,z) and the τ_(xy,yz,xz) and determining if any residuals exist, the rock displacement parameters u, v, and w representing accurate rock displacement parameters for said reservoir when the residuals are substantially equal to zero; and (e) determining said subsidence in said reservoir from said accurate rock displacement parameters determined from step (d), by solving a final set of failure criterion equations, which comprise the residuals and any derivatives, to determine a set of rock plastic displacements, said rock plastic displacements forming a part of said rock displacement parameters u, v, and w.
 8. The program storage device of claim 7, wherein the determining step (b) comprises the step of: determining said ε_(x,y,z) (the ‘x,y,z elongation strains’) and said γ_(xy,yz,zx) (the ‘shear strains’) from said rock displacement parameters u, v, and w by using the following equations: $\begin{matrix} {ɛ_{x} = \frac{\partial u}{\partial x}} \\ {ɛ_{y} = \frac{\partial v}{\partial y}} \\ {ɛ_{z} = \frac{\partial w}{\partial z}} \\ {\gamma_{xy} = {\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}}} \\ {\gamma_{yz} = {\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}}} \\ {\gamma_{zx} = {\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}}} \end{matrix}$


9. The program storage device of claim 8, wherein the determining step (c) comprises the step of: determining said σ_(x,y,z) (the ‘elastic normal rock stress in x,y,z directions’) and said τ_(xy,yz,xz) (the ‘elastic shear stress’) from said ε_(x,y,z) and γ_(xy,yz,zx) by using the following equations: $\begin{matrix} {\sigma_{x} = {{2G\quad ɛ_{x}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\ {\sigma_{y} = {{2G\quad ɛ_{y}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\ {\sigma_{z} = {{2G\quad ɛ_{z}} + {\lambda \left( {ɛ_{x} + ɛ_{y} + ɛ_{z}} \right)}}} \\ {\tau_{xy} = {G\quad \gamma_{xy}}} \\ {\tau_{yz} = {G\quad \gamma_{yz}}} \\ {\tau_{zx} = G_{\gamma_{zx}}} \end{matrix}$


10. The program storage device of claim 9, wherein the solving step (d) comprises the steps of: solving said set of rock momentum balance differential equations from said σ_(x,y,z) and said τ_(xy,yz,xz), wherein said rock momentum balance differential equations have a left-hand side on a left side of an equal sign and a right-hand side on a right side of said equal sign, said rock momentum balance differential equations including, $\begin{matrix} {{\frac{\partial\sigma_{x}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} + F_{x} + P_{x}} = 0} \\ {{\frac{\partial\sigma_{y}}{\partial y} + \frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{zy}}{\partial z} + F_{y} + P_{y}} = 0} \\ {{\frac{\partial\sigma_{z}}{\partial z} + \frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + F_{z} + P_{z}} = 0} \end{matrix}$

said rock displacement parameters u, v, and w estimated in step (a) representing said accurate rock displacement parameters for said reservoir when a left-hand side of said rock momentum balance differential equations substantially equal a right hand side of said rock momentum balance differential equations,
 11. The program storage device of claim 10, wherein the solving step (e) comprises the steps of using the said σ_(x,y,z) and said τ_(xy,yz,xz) to determine whether a failure criterion has been exceeded in one of a set of staggered grids.
 12. The program storage device of claim 11, wherein, when said failure criterion has been exceeded in said one of said set of staggered grids, the solving step (e) further comprises the step of: setting up a final set of equations, each equation being a failure criterion approximately equal to zero, and determining a set of parameters up, vp, wp which are plastic displacements resulting from a failure of the rock, said plastic displacements forming a part of said subsidence and being determined from said accurate rock displacement parameters for said reservoir. 